#!/usr/bin/perl -w
use strict;
use FindBin;
use Getopt::Long;
use lib "/net/cpp-group/Leo/bin";
use parse_fasta;
use Cwd;

#88888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Usage

#88888888888888888888888888888888888888888888888888888888888888888888888888888888
my $usage = <<"USAGE";

USAGE:

    site_specific_pipeline_slr.pl --cdna cdna.fa --prot prot.seq --f3x4
                            [ --msf --clustal --regex "(.*)" --verbose --help]

        This script accepts protein sequences in either CLUSTAL or MSF or
        FASTA format (default) and cDNA in FASTA format.

        The matching IDs from the protein and cDNA sequence files can be
        identified using a regular expression. The default is to take the
        first continuous stretch of non-space characters

		Results in results.parameters

	Check:
		Always check for errors in the pipeline (All files ending in '.err')
		   ***ESPECIALLY*** 
		   mapping.err
		   
	Happy Hunting.
				  .... Leo (Goodstadt) 1/11/2004 to 27/05/2005

USAGE

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Get options

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
# mandatory


# optional parameters
my $help = undef;
my $regex = '(\\S+)';

my $prot_file;
my $cdna_file;
my $msf;
my $f3x4;
my $clustal;
my $verbose;
my $vverbose;
#{
#    open CMDLINE, ">cmd.line" or die "Error:\n\tCould not open cmd.line\n";
#    print CMDLINE join (" ", @ARGV), "\n";
#    close CMDLINE;
#}
GetOptions(
			'prot=s'	=> \$prot_file,
			'cdna=s'	=> \$cdna_file,
			'regex=s'	=> \$regex,
			'f3x4'		=> \$f3x4,
			'msf'		=> \$msf,
			'clustal'	=> \$clustal,
			'verbose'	=> \$verbose,
			'vverbose'	=> \$vverbose,
			'help'		=> \$help,
			);

die $usage if ($help);
die $usage unless ($prot_file && $cdna_file);
die "Error\n\t:[$prot_file] does not exist.\n" unless (-e $prot_file);
die "Error\n\t:[$cdna_file] does not exist.\n" unless (-e $cdna_file);
my $verbose_cmd = $vverbose ? "--verbose" :"";
unless ($regex =~ /^".*"$/)
{
	$regex = '"'.$regex.'"';
}

$f3x4 = defined ($f3x4) ? 1 : 0;

$verbose = 1 if ($vverbose);
my $bin_dir = $FindBin::Bin;
$bin_dir = substr($bin_dir, 0, -1) if substr($bin_dir, -1, 1) eq '/';

sub my_mkdir($)
{
	my ($dir_name) = @_;
	if (!-e $dir_name)
	{
		mkdir $dir_name
				or die "Error\n\tCould not create the directory [$dir_name]\n$!\n";
	}
	elsif (! -d  "$dir_name")
	{
		die "Error\n\t[$dir_name] is a file not a directory.\n$!\n";
	}

}
sub run_cmd
{
	#remove beginning white space and new lines
	my $cmd = $_[0];
	$cmd =~tr/\t\n/ /d;
	$cmd =~ s/\s+/ /gm;
	if ($_[1])
	{
		print
#STDERR
 "\tat ", getcwd(), "\n";
		print
#STDERR
 "\t", $cmd, "\n";
	}
	system $cmd and die "Error:\n\tFailed [ $cmd ]\n\t$?\n\t$!\n";
}

sub print_site_specific_ctl;
sub print_pairwise_ctl;
print "\n\n" if $verbose;

#
#	Prepare FASTA format protein sequence file
#
if ($msf)
{
	print
#STDERR
 "MSF --> FASTA ...\n" if $verbose;
run_cmd(<<"CMD");
	$bin_dir/msf_to_fasta.pl $prot_file
		> $prot_file.fa
CMD
}
elsif ($clustal)
{
	print
#STDERR
 "CLUSTAL --> FASTA ...\n" if $verbose;
run_cmd(<<"CMD");
	$bin_dir/clustal_to_fasta.pl $prot_file
		> $prot_file.fa
CMD
}
else
{
run_cmd(<<"CMD");
	cp $prot_file $prot_file.fa
CMD
}

#
#	Map cDNA onto aa
#
{
	print
#STDERR
 "Using protein alignment to align corresponding cDNA ...\n" if $verbose;
run_cmd(<<"CMD");
	$bin_dir/map_cdna_onto_aa
		--peptide_sequences $prot_file.fa
        --cdna_sequences $cdna_file
        --peptide_output $prot_file.mapped.fa
        --cdna_output $cdna_file.mapped.fa
        --regex $regex
        --err_log mapping.err
		--x_for_aa_mismatches X
		$verbose_cmd
CMD
	unlink ("$prot_file.fa");
}

#
#	Prepare cDNA onto aa
#
{
	print
#STDERR
 "Convert mapped cDNA to PAML format ...\n" if $verbose;
run_cmd(<<"CMD");
	$bin_dir/seq_extract
		--wrap_line 0
		--output_format p
		$verbose_cmd
		--e extract.err
		< $cdna_file.mapped.fa
		> seq.tmp
CMD
}

#
#  -------Create pairwise Ks directory
#
#
#  -------make directories
#
{
	print
#STDERR
 "Make directories for codeml analyses ...\n" if $verbose;
	my_mkdir("codeml.pairwise");
	my_mkdir ("slr");
}



#
#  -------Create PAML control files ..
#
print
#STDERR
 "Create PAML control files ...\n" if $verbose;
print_site_specific_ctl();
print_pairwise_ctl();





#
#  -------Calculate pairwise Ks to build tree
#
#=pod
{
	print
#STDERR
 "Calculate pairwise Ks to build tree ...\n" if $verbose;
run_cmd(<<"CMD");
	$bin_dir/seq_pairs_kaks
		-aA
		-t codeml.pairwise
		-c $cdna_file.mapped.fa
		--err_log kaks.pairwise.err
		-p paml.err
		$verbose_cmd
		> kaks.results
CMD

print "Filter out Ks for build tree ...\n" if $verbose;
run_cmd(<<"CMD");
	$bin_dir/seq_pairs_filter_kaks_results
		-S
		-R filter.summary
		-r 0
		-f 0
		-s 0
		-t 0
		-n 0
		-l 0
		-e filter.err
		$verbose_cmd
		< kaks.results
		> codeml.pairwise/pairwise.dS
CMD


#	chdir ("codeml.pairwise");
#	my $cmd = "$bin_dir/../../tools/paml/codeml pairwise_kaks.ctl |";
#	my $cmd = "$bin_dir/../../tools/paml/codeml pairwise_kaks.ctl |";
#	my $cmd = "$bin_dir/unbuffer $bin_dir/../../tools/paml/codeml pairwise_kaks.ctl |";
#	open(RUN1, $cmd)
#		or die "Error\n\tCould not run [$cmd]\n$!\n";
#
#	while (<RUN1>)
#	{
#		/^\s+(.*?)\s+$/;
#		print
##STDERR
# "\t$1\n" if ($1 && length($1) && $verbose);
#	}
#	chdir "..";
#	unless (-e "codeml.pairwise/2ML.dS")
#	{
#		die "Error:\n\tPairwise codeml failed to run successfully.\n".
#				"The file are in the [codeml.pairwise] sub-directory.\n";
#	}
#
#	# rename to pairwise dS before feeding into fitch
#	run_cmd(<<"CMD");
#	mv codeml.pairwise/2ML.dS codeml.pairwise/pairwise.dS
#CMD
}


#
#  -------Build fitch tree from ks values
#
{
	print
#STDERR
 "Build fitch tree from ks values ...\n" if $verbose;

	# run fitch
	run_cmd(<<"CMD");
	$bin_dir/tree_fitch -i p
			$verbose_cmd
			< codeml.pairwise/pairwise.dS
			> fitch.tree
CMD
}


#
#  -------site specific analyses
#
{

	# 0.03
	print "Site specific analyses with slr ...\n" if $verbose;
	chdir("slr");
	system("$bin_dir/tree_to_slr_format.pl < ../fitch.tree --sequence ../seq.tmp > slr.tree");
	system("cp ../seq.tmp slr.seq");
	print STDERR "run Slr:\n";
	print STDERR "$bin_dir/tools/slr/bin/Slr slr.ctl > output.slr\n";
	system("$bin_dir/tools/slr/bin/Slr slr.ctl > output.slr");
	system("mv slr.results ..");
	system("mv output.slr ../slr.output");
	chdir("..");
}

#=cut

{
	print
#STDERR
 "Analyse results ...\n" if $verbose;

	# summarise SLR output
	{
		open (SLR_STDOUT, "slr.output") or die;
		open (SUMMARY, ">output.summary") or die;
		while (<SLR_STDOUT>)
		{
			last if /^\*\*\*/;
		}
		
		$_ = <SLR_STDOUT>; # ignore next

		# tree
		{
			$_ = <SLR_STDOUT>;
			chomp;
			print SUMMARY "$_\n";
		}
		
		# lnL
		{
			my $lnL = <SLR_STDOUT>;
			$lnL =~ /# lnL = ([0-9\.\+\-\e]+).*/ or die;
			print SUMMARY "$1\tlnL\n";
		}

		# kappa / omega
		{
			my $line = <SLR_STDOUT>;
			$line =~ /Kappa = ([0-9\.\+\-\e]+) Omega = ([0-9\.\+\-\e]+)/ or die;
			print SUMMARY "$1\tkappa\n";
			print SUMMARY "$2\tomega\n";
		}

		# tree parameters
		{
			my $tree = <SLR_STDOUT>;
			$tree =~ /Tree length = ([0-9\.\+\-\e]+), average branch length = ([0-9\.\+\-\e]+) \(min=([0-9\.\+\-\e]+), max=([0-9\.\+\-\e]+)\)/ or die;
			print SUMMARY
			"$1\ttree length\n",
			"$2\taverage branch length\n",
			"$3\tminimum branch length\n",
			"$4\tmaximum branch length\n";
		}
	}

	run_cmd(<<"CMD");
	$bin_dir/parse_slr_results.pl
			--significant_residues  results.significant.fa
			$verbose_cmd
			slr.results
			>> output.summary
CMD

	
	system("more output.summary");

	
}



print "All finished...\n\n";








#
#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	parameter files
#







sub print_pairwise_ctl
{
	open PAIR_CTL, ">codeml.pairwise/pairwise_kaks.ctl"
			or die "Error\n\t:Could not write to pairwise_kaks.ctl\n$!\n";
	print PAIR_CTL <<"CMD";
seqfile = ../seq.tmp
outfile = results
icode = 0	  			# universal code
noisy = 0
verbose = 0
runmode = -2  			# pairwise
seqtype = 1	  			# codons
CodonFreq = 2 			# F3X4; Could be 3 for codon table
aaDist = 0	  			# ??? equal Could be geometric or linear or G1974 Miyata
model = 0	  			# one dN/dS over tree
NSsites = 0	  			# one dN/dS over sequence
Mgene = 0		    	# ??? rates; 1: separate
fix_omega = 0
omega = .4
fix_alpha = 1		 	# gamma distribution shape set to 1 for single rate
alpha = 0		     	# gamma distribution shape set to 1 for single rate
fix_rho = 1
rho = 0
Malpha = 0		     	# ??? different alphas for genes
ncatG = 1		     	# # of categories in dG of NSsites models5D
clock = 0		     	# ??? no clock 1:clock
getSE = 1
RateAncestor = 1	 	# ??? rates alpha > 0 or ancestral rates 1 or 2
Small_Diff = .5e-6
method = 0		     	# simultaneous, 1: one branch at a time
fix_blength = 0		 	# ignore branch lengths from supplied ree
fix_kappa = 0
kappa = 2
CMD
}

sub print_site_specific_ctl()
{
	my ($start_omega) = @_;
	my $file = ">slr/slr.ctl";
	open SLR_CTL, $file
			or die "Error\n\t:Could not write to $file\n$!\n";
print SLR_CTL <<"CMD";
seqfile: slr.seq
treefile: slr.tree
outfile: slr.results
positive_only: 0
reoptimize: 1
codonf: $f3x4
CMD
}

print_site_specific_ctl();

sub print_fitch_parameters
{
	my $file = ">fitch.parameters";
	open FITCH_CTL, $file
			or die "Error\n\t:Could not write to $file\n$!\n";
	print FITCH_CTL <<"CMD";
build_tree.dS
0
7
L
J
33
50
2
Y
CMD
	close FITCH_CTL;
}



